# Specify the columns to check for NA values# List of columns to check for NA valuescolumns_to_check <-c("DroughtTrt", "DroughNet_plotID", "ageClass", "siteID", "species", "plant_height", "SLA", "mean_thickness", "dry_mass_g_original", "LDMC", "leaf_area", "wet_mass_g", "dry_mass_g")# Check for rows with any NA values in the specified columnsna_rows <-apply(droughtnet[columns_to_check], 1, function(row) any(is.na(row)))# Create a data frame showing how many NA values there are in each rowna_row_summary <- droughtnet[na_rows, columns_to_check]na_row_summary$row_na_count <-rowSums(is.na(na_row_summary))# Summarize the number of rows by the count of NAs they containna_row_summary <- na_row_summary %>%group_by(row_na_count) %>%summarise(n =n(), .groups ='drop')# Create a summary showing the co-occurrence of NA values across the specified columnsna_co_occurrence <- droughtnet %>%summarise(across(all_of(columns_to_check), ~sum(is.na(.)))) %>%pivot_longer(cols =everything(), names_to ="variable", values_to ="na_count") %>%arrange(desc(na_count))# Print the summariesprint(paste("Total rows with any NA values across checked columns:", sum(na_rows)))
[1] "Total rows with any NA values across checked columns: 50"
# Create the new data frame with selected columnsdroughtnet_data2 <-select(droughtnet, "envelope_ID", "siteID", "species", "ageClass", "DroughtTrt", "DroughNet_plotID", "leaf_age", "calluna_shoot_type", "plant_height", "mean_thickness", "SLA", "LDMC", "plant_nr")view(droughtnet_data2)# Specify the columns to check for NA valuescolumns_to_check <-c("plant_height", "mean_thickness", "SLA", "LDMC")# Remove rows with any NA values in the specified columnsdroughtnet_data2_clean <- droughtnet_data2[!rowSums(is.na(droughtnet_data2[columns_to_check])), ]view(droughtnet_data2_clean)# Calculate the total number of NAs removed (for the specified columns only)# This calculation considers rows removed rather than individual NAstotal_NAs_removed <-nrow(droughtnet_data2) -nrow(droughtnet_data2_clean)# Calculate the total number of entries remainingtotal_entries_remaining <-nrow(droughtnet_data2_clean)# Print out the total NAs removed and total entries remainingcat("Total NAs removed:", total_NAs_removed, "\n")
Analysis of young leaves for the four focal species
# Specify the species of interest, including Calluna vulgarisspecies_of_interest <-c("Empetrum nigrum", "Vaccinium myrtillus", "Vaccinium vitis-idaea", "Calluna vulgaris")# Filter the data to include only the specified species where leaf_age is 'young',# and for Calluna vulgaris, additionally check that calluna_shoot_type is 'short'filtered_data <- droughtnet_data2_clean %>%filter(species %in% species_of_interest, leaf_age =="young") %>%filter(# Include Calluna vulgaris from Tjøtta regardless of calluna_shoot_type (species =="Calluna vulgaris"& siteID =="Tjøtta") |# Include Calluna vulgaris from Lygra only if calluna_shoot_type is Short (species =="Calluna vulgaris"& siteID =="Lygra"& calluna_shoot_type =="Short") |# Include other species of interest from both sites (species !="Calluna vulgaris") )view(filtered_data)
D. plots for young leaves for the traits
filtered_data <- filtered_data %>%filter(siteID %in%c("Tjøtta", "Lygra"), DroughtTrt %in%c("Amb (0)", "Ext (90)"), ageClass %in%c("Pioneer", "Mature"))ggplot(filtered_data, aes(x = species, y = plant_height, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +# Removed leaf_age from facet_gridtheme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="Plant Height vs Treatment by Species for Young Leaves",x ="Species",y ="Plant Height") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey")) # Assigning colors to treatments
ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, vjust =0.5, hjust =1), # vertical adjustment for x axis labelsplot.caption =element_text(hjust =0.5)) +# Center the captionlabs(title ="SLA vs Treatment by Species",x ="Species",y ="Log (SLA)",fill ="Drought Treatment",caption ="Figure 2: SLA Variation in Dwarf Shrubs Under Drought Conditions Across Different Ages and Sites") +scale_fill_manual(values =c("Amb (0)"="orange", "Ext (90)"="darkgrey")) # Assigning colors to treatments
ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, vjust =0.5, hjust =1), # vertical adjustment for x axis labelsplot.caption =element_text(hjust =0.5)) +# Center the captionlabs(title ="LDMC vs Treatment by Species",x ="Species",y ="LDMC",fill ="Drought Treatment",caption ="Figure 3: LDMC Variation in dwarf Shrubs Under drought conditions across different successional phases and sites") +scale_fill_manual(values =c("Amb (0)"="darkblue", "Ext (90)"="lightblue")) # Assigning colors to treatments
filtered_data <- filtered_data %>%filter(siteID %in%c("Tjøtta", "Lygra"), DroughtTrt %in%c("Amb (0)", "Ext (90)"), ageClass %in%c("Pioneer", "Mature"))ggplot(filtered_data, aes(x = species, y = mean_thickness, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +# Removed leaf_age from facet_gridtheme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="Plant Height vs Treatment by Species for Young Leaves",x ="Species",y ="mean_thickness") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey")) # Assigning colors to treatments
Trait 1- LDMC- models
# fitting a linear mixed-effects model using lmer, plant nr is the nested random effect# Ensure variables are treated as categorical if they represent categoriesfiltered_data$DroughNet_plotID <-as.factor(filtered_data$DroughNet_plotID)# Model with all the main effects for LDMCmodel_main_effects_LDMC <-lmer(LDMC ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects_LDMC)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
LDMC ~ species + DroughtTrt + ageClass + siteID + (1 | DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: 8228.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.7757 -0.5602 0.0437 0.5846 5.7000
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 300.3 17.33
DroughNet_plotID (Intercept) 171.6 13.10
Residual 2643.3 51.41
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 349.198 8.257 26.279 42.292 < 2e-16 ***
speciesEmpetrum nigrum 8.147 5.276 721.162 1.544 0.122988
speciesVaccinium myrtillus 2.088 5.481 700.600 0.381 0.703352
speciesVaccinium vitis-idaea -21.253 5.136 690.820 -4.138 3.93e-05 ***
DroughtTrtExt (90) -2.987 7.649 19.272 -0.390 0.700462
ageClassPioneer 30.674 7.659 19.413 4.005 0.000731 ***
siteIDTjøtta 4.638 7.649 19.320 0.606 0.551346
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
# The intercept for LDMC, representing calluna at reference level, is significantly high at 359.3762, indicating a strong baseline effect on LDMC.# Species Impact: The species Vaccinium vitis-idaea significantly reduces LDMC by 30.0715 compared to the baseline species- unique biological characteristics.# Drought Treatment: Drought treatment does not significantly affect LDMC, as indicated by a p-value of 0.160270, suggesting drought stress may not be a critical factor for LDMC variation.# Age Class Influence: Plants classified under the pioneer age class show a significant increase in LDMC (26.5932), pointing towards age or developmental stage as a determinant of LDMC.# Site Variation: Tjøtta, does not significantly influence LDMC, with a p-value of 0.152627, indicating that location-based environmental factors may not play a significant role in LDMC variationperformance::check_model(model_main_effects_LDMC, detrend =FALSE)
#The Posterior Predictive Check looks reasonable, as the observed data seem to match the model-predicted data well. #The Linearity plot suggests a linear relationship # The Homogeneity of Variance plot does not show an obvious fan shape which would indicate heteroscedasticity, so it seems the assumption of homoscedasticity is met. # The Influential Observations plot indicates a few points with high leverage or large residuals,-but not worrisome as they are extreme values # The Collinearity plot suggests no issues with collinearity since all VIF values are well below the common threshold. # The Normality of the Residuals plot indicates some deviation from normality, especially in the tails. # The Normality of Random Effects plots show some slight deviations from normality #overall, the model fits the data well
fitting model to the ggplot
# Get fitted values directly from the modelfitted_values <-predict(model_main_effects_LDMC, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
interaction models
#species * DroughtTrtmodel_species_droughtTrt <-lmer(LDMC ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model_species_droughtTrt)
# intercept(calluna) LDMC at 365.795 indicate a significant base level of LDMC before considering species variation or drought effects.# Species-Specific Effects:# Empetrum nigrum notably increases LDMC by 23.225 units compared to calluna demonstrating- could be a positive adaptation # Vaccinium myrtillus also shows a positive but not statistically significant effect on LDMC, suggesting a mild or variable adaptation compared to calluna# Vaccinium vitis-idaea significantly reduces LDMC by 29.776 units, highlighting its distinct physiological or ecological traits that might contribute to lower LDMC# Drought Treatment Response:# drought trt does not significantly alter LDMC on its own# Interaction with Drought:# Both Empetrum nigrum and Vaccinium myrtillus exhibit significant negative interactions with drought treatment, indicating a decrease in LDMC under drought condition- suggests these species could be sensitive to drought, affecting their LDMC.# Vaccinium vitis-idaea shows no significant interaction with drought#validate modelperformance::check_model(model_species_droughtTrt, detrend =FALSE)
Fitting model to plot
# Get fitted values directly from the modelfitted_values <-predict(model_species_droughtTrt, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_values# Now create your plotplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
# species * age Class (h2)#this will help adress the hypothesis by examining weather the effect of drought on trait varies between age classes#significant interaction supports the hypothesis that mat and pio respond differently to droght conditions, potentially more consrvative traits in mat than pioLDMC_species_ageClass <-lmer(LDMC ~ species * ageClass + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(LDMC_species_ageClass)
# Get fitted values directly from the modelfitted_values <-predict(LDMC_species_ageClass, re.form =NA) filtered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
# Drought Treatment * siteID (related to h3)#will help to show impact between north and south# significant interaction would indicate regional variability in how drought influences trait- hence supportinng the hypothesis, otherwise then..model1_species_siteID <-lmer(LDMC ~ species * siteID + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model1_species_siteID)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: LDMC ~ species * siteID + (1 | DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: 8150.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.6641 -0.5989 0.0624 0.5847 6.2545
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 361.0 19.00
DroughNet_plotID (Intercept) 326.1 18.06
Residual 2347.7 48.45
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value
(Intercept) 369.128 7.766 43.839 47.529
speciesEmpetrum nigrum 25.754 7.581 746.502 3.397
speciesVaccinium myrtillus 4.138 7.650 693.341 0.541
speciesVaccinium vitis-idaea -56.769 7.087 687.616 -8.011
siteIDTjøtta -7.349 10.908 42.636 -0.674
speciesEmpetrum nigrum:siteIDTjøtta -26.068 10.087 725.243 -2.584
speciesVaccinium myrtillus:siteIDTjøtta -1.468 10.402 691.426 -0.141
speciesVaccinium vitis-idaea:siteIDTjøtta 64.993 9.721 685.887 6.686
Pr(>|t|)
(Intercept) < 2e-16 ***
speciesEmpetrum nigrum 0.000717 ***
speciesVaccinium myrtillus 0.588777
speciesVaccinium vitis-idaea 4.89e-15 ***
siteIDTjøtta 0.504115
speciesEmpetrum nigrum:siteIDTjøtta 0.009952 **
speciesVaccinium myrtillus:siteIDTjøtta 0.887791
speciesVaccinium vitis-idaea:siteIDTjøtta 4.77e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
# Get fitted values directly from the modelfitted_values <-predict(model1_species_siteID, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values plot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 1: Shows LDMC of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
#model with all the main effectsmodel_main_effects <-lmer(plant_height ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: plant_height ~ species + DroughtTrt + ageClass + siteID + (1 |
DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: 4558.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.7321 -0.5753 -0.0484 0.6013 4.0400
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 3.410 1.847
DroughNet_plotID (Intercept) 6.381 2.526
Residual 19.915 4.463
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 26.2221 1.1926 22.1651 21.987 < 2e-16 ***
speciesEmpetrum nigrum -5.5906 0.4611 716.9459 -12.123 < 2e-16 ***
speciesVaccinium myrtillus -9.6354 0.4782 696.0138 -20.150 < 2e-16 ***
speciesVaccinium vitis-idaea -10.8049 0.4470 693.9222 -24.174 < 2e-16 ***
DroughtTrtExt (90) 3.3234 1.1615 19.9300 2.861 0.009673 **
ageClassPioneer -4.3616 1.1619 19.9685 -3.754 0.001253 **
siteIDTjøtta -5.1761 1.1614 19.9289 -4.457 0.000244 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
#the model does not show any normality of residuals and also there is heterogeneity of variance#need to modify by transforming the fitted values#plant heightneeds to be log transformed to achieve homoscedasticity and linearitymodel_main_effects3 <-lmer(log(plant_height) ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(plant_height) ~ species + DroughtTrt + ageClass + siteID +
(1 | DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: 188.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.97832 -0.62162 0.01327 0.64119 2.77955
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 0.01764 0.1328
DroughNet_plotID (Intercept) 0.03083 0.1756
Residual 0.06010 0.2452
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.24672 0.08126 21.26656 39.956 < 2e-16
speciesEmpetrum nigrum -0.27831 0.02548 713.70008 -10.921 < 2e-16
speciesVaccinium myrtillus -0.53579 0.02635 696.22629 -20.331 < 2e-16
speciesVaccinium vitis-idaea -0.60639 0.02463 695.59109 -24.624 < 2e-16
DroughtTrtExt (90) 0.19009 0.07985 19.82276 2.381 0.02744
ageClassPioneer -0.26757 0.07987 19.85403 -3.350 0.00321
siteIDTjøtta -0.28991 0.07985 19.82382 -3.631 0.00168
(Intercept) ***
speciesEmpetrum nigrum ***
speciesVaccinium myrtillus ***
speciesVaccinium vitis-idaea ***
DroughtTrtExt (90) *
ageClassPioneer **
siteIDTjøtta **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
# Species: Empetrum nigrum, Vaccinium myrtillus, and Vaccinium vitis-idaea all show significant decreases in height compared to Calluna, indicating notable species-specific differences.# Drought Treatment: drought positively affects height, suggesting an increase under drought conditions for all species# Age Class: The pioneer age class is associated with a significant decrease in height, indicating lower levels compared to more mature.# Site: The Tjøtta site is linked to a decrease in height, highlighting the influence of location on height#validate the model- check for linearity, homoscedasticity, outliers by plotting the resudualsperformance::check_model(model_main_effects3, detrend =FALSE)
The diagnostic plots indicate that the model has a generally good fit. The predictions match the observed data well, the residuals show linearity and homogeneity of variance, though slight deviation, and they are also normally distributed. The random effects appear to be normally distributed
fitt model to plot
# Get fitted values directly from the modelfitted_values <-predict(model_main_effects3, re.form =NA) # 'NA' to exclude random effects# Exponentiate the fitted values to transform back to the original scalefiltered_data$fitted <-exp(fitted_values)plot <-ggplot(filtered_data, aes(x = species, y = plant_height, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="Plant Height vs Treatment by Species",x ="Species", y ="Plant Height", # Updated y-axis labelcaption ="Figure 1: Shows plant height of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer and mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)
interaction models
#models for two way interactions#species * drought trtmodel_species_drought2 <-lmer(log(plant_height) ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_species_drought2)
# Compared to Calluna, the other species significantly reduce the height, with Vaccinium species showing the largest decreases.# Drought treatment on its own does not significantly affect theight, but its interaction with species (Vaccinium myrtillus and Vaccinium vitis-idaea) leads to significant increases, suggesting these species may respond positively to drought conditions in terms of thier height#plot diagnostic plotsperformance::check_model(model_species_drought2, detrend =FALSE)
# Get fitted values directly from the modelfitted_values <-predict(model_species_drought2, re.form =NA) # Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y =log(plant_height), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="red", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="plant height vs Treatment by Species",x ="Species", y ="LDMC") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
# Get fitted values directly from the modelfitted_values <-predict(model_height_species_ageclass, re.form =NA) # Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_values# Now create your plotplot <-ggplot(filtered_data, aes(x = species, y =log(plant_height), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="red", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="plant height vs Treatment by Species",x ="Species", y ="plant height") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
# model for sla with all the fixed effectsmodel_main_effects <-lmer(SLA ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
SLA ~ species + DroughtTrt + ageClass + siteID + (1 | DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: 7482.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.5193 -0.5305 -0.0876 0.3306 6.3637
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 104.28 10.212
DroughNet_plotID (Intercept) 39.66 6.298
Residual 996.03 31.560
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 87.651 4.588 29.915 19.105 < 2e-16 ***
speciesEmpetrum nigrum 42.415 3.233 730.780 13.121 < 2e-16 ***
speciesVaccinium myrtillus 155.145 3.359 715.032 46.183 < 2e-16 ***
speciesVaccinium vitis-idaea 45.801 3.151 706.068 14.537 < 2e-16 ***
DroughtTrtExt (90) -3.031 4.171 20.348 -0.727 0.475771
ageClassPioneer -10.719 4.179 20.569 -2.565 0.018221 *
siteIDTjøtta -16.294 4.172 20.439 -3.905 0.000849 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
##the model does not fit the data well#need to transform the observed values to achive normality and homoogeneity of variance#log transform SLA to achive normality and homodeneity of variancesmodel_main_effects4 <-lmer(log(SLA) ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects4)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
log(SLA) ~ species + DroughtTrt + ageClass + siteID + (1 | DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: -133.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.2558 -0.5792 -0.0635 0.5222 3.2135
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 0.005324 0.07297
DroughNet_plotID (Intercept) 0.001795 0.04236
Residual 0.042768 0.20680
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.373178 0.031036 28.010849 140.907 < 2e-16
speciesEmpetrum nigrum 0.457574 0.021221 723.224265 21.563 < 2e-16
speciesVaccinium myrtillus 1.134496 0.022039 703.719975 51.476 < 2e-16
speciesVaccinium vitis-idaea 0.455653 0.020664 694.265373 22.050 < 2e-16
DroughtTrtExt (90) -0.003061 0.028388 19.502029 -0.108 0.915241
ageClassPioneer -0.081935 0.028439 19.716246 -2.881 0.009325
siteIDTjøtta -0.119785 0.028394 19.589008 -4.219 0.000439
(Intercept) ***
speciesEmpetrum nigrum ***
speciesVaccinium myrtillus ***
speciesVaccinium vitis-idaea ***
DroughtTrtExt (90)
ageClassPioneer **
siteIDTjøtta ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
# Calluna exhibits a significant baseline SLA value of 4.327885, serving as the reference for other species comparisons.# Species Effects:# Empetrum nigrum Significantly increases SLA by 0.519248 compared to Calluna.# Vaccinium myrtillus Shows a substantial increase in SLA by 1.195279, the largest among the species.# Vaccinium vitis-idaea also significantly increases SLA by 0.513608, similar to Empetrum nigrum.# All species effects are highly significant, indicating distinct SLA characteristics compared to Calluna.# Drought treatment effect on SLA is negligible (-0.001773) and not statistically significant, suggesting drought has minimal/no impact on SLA across the specie# Age Class (Pioneer)indicates a significant decrease in SLA (-0.083333), suggesting pioneer plants have a lower SLA compared to more mature stages.# Site (Tjøtta): Associated with a significant decrease in SLA (-0.152224), indicating environmental or locational factors at Tjøtta contribute to lower SLA values#validate the model- check for linearity, homoscedasticity, outliers by plotting the resudualsperformance::check_model(model_main_effects4, detrend =FALSE)
Overall, the model shows a generally good fit with some concerns:
The posterior predictive check indicates a good fit. The linearity assumption seems to be met. There is a potential issue with homoscedasticity, as the Homogeneity of Variance plot suggests that variance slightly increases with fitted values. Influential observations do not appear to be a problem. The residuals and random effects are approximately normally distributed
fit model onto the plot
# Get fitted values directly from the modelfitted_values <-predict(model_main_effects4, re.form =NA) # Exponentiate the fitted values to transform back to the original scalefiltered_data$fitted <-exp(fitted_values)# Now create your plot with a caption at the bottomplot <-ggplot(filtered_data, aes(x = species, y = SLA, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +# Use exponentiated fitted valuesfacet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA", caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer and mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
Interaction models for SLA
#species * droughtTrtsla_species_droughtTrt <-lmer(log(SLA) ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(sla_species_droughtTrt)
In summary, the model seems to have a good fit overall. high VIF for the ‘species:DroughtTrt’ interaction term is acceptable, and slight deviations from normality in the residuals. As for the residuals, the minor deviations from normality .
# Get fitted values directly from the modelfitted_values <-predict(sla_species_droughtTrt, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)
#species * ageclasssla_species_ageClass <-lmer(log(SLA) ~ species * ageClass + (1| DroughNet_plotID), data = filtered_data)summary(sla_species_ageClass)
# Get fitted values directly from the modelfitted_values <-predict(sla_species_ageClass, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <-exp(fitted_values) # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y = SLA, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)
# species * siteIDsla_species_siteID <-lmer(log(SLA) ~ species * siteID + (1| DroughNet_plotID), data = filtered_data)summary(sla_species_siteID)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(SLA) ~ species * siteID + (1 | DroughNet_plotID)
Data: filtered_data
REML criterion at convergence: -153.7
Scaled residuals:
Min 1Q Median 3Q Max
-4.8127 -0.5784 -0.0517 0.5263 3.3217
Random effects:
Groups Name Variance Std.Dev.
DroughNet_plotID (Intercept) 0.004578 0.06766
Residual 0.043666 0.20896
Number of obs: 765, groups: DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df
(Intercept) 4.30581 0.02861 58.33859
speciesEmpetrum nigrum 0.43005 0.03168 754.09060
speciesVaccinium myrtillus 1.16166 0.03258 751.05284
speciesVaccinium vitis-idaea 0.57441 0.03024 739.17435
siteIDTjøtta -0.05570 0.04006 56.13847
speciesEmpetrum nigrum:siteIDTjøtta 0.02295 0.04258 747.84125
speciesVaccinium myrtillus:siteIDTjøtta -0.06125 0.04436 747.92330
speciesVaccinium vitis-idaea:siteIDTjøtta -0.23071 0.04152 737.49954
t value Pr(>|t|)
(Intercept) 150.506 < 2e-16 ***
speciesEmpetrum nigrum 13.575 < 2e-16 ***
speciesVaccinium myrtillus 35.653 < 2e-16 ***
speciesVaccinium vitis-idaea 18.998 < 2e-16 ***
siteIDTjøtta -1.390 0.170
speciesEmpetrum nigrum:siteIDTjøtta 0.539 0.590
speciesVaccinium myrtillus:siteIDTjøtta -1.381 0.168
speciesVaccinium vitis-idaea:siteIDTjøtta -5.557 3.83e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
# Get fitted values directly from the modelfitted_values <-predict(sla_species_siteID, re.form =NA) # 'NA' to exclude random effectsfiltered_data$fitted <- fitted_values # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)
#general model#model with all the main effectsmodel_main_mean_thickness_effects <-lmer(mean_thickness ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data) summary(model_main_mean_thickness_effects )
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: mean_thickness ~ species + DroughtTrt + ageClass + siteID + (1 |
DroughNet_plotID/plant_nr)
Data: filtered_data
REML criterion at convergence: -1873.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.1179 -0.5755 -0.1043 0.3895 4.9624
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 0.0001841 0.01357
DroughNet_plotID (Intercept) 0.0001764 0.01328
Residual 0.0044818 0.06695
Number of obs: 765, groups:
plant_nr:DroughNet_plotID, 82; DroughNet_plotID, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.837e-01 8.906e-03 3.206e+01 43.087 < 2e-16
speciesEmpetrum nigrum -2.059e-01 6.807e-03 7.302e+02 -30.250 < 2e-16
speciesVaccinium myrtillus -2.862e-01 7.091e-03 7.146e+02 -40.358 < 2e-16
speciesVaccinium vitis-idaea -1.169e-01 6.659e-03 6.995e+02 -17.558 < 2e-16
DroughtTrtExt (90) -7.296e-05 7.949e-03 2.021e+01 -0.009 0.993
ageClassPioneer 7.590e-04 7.961e-03 2.036e+01 0.095 0.925
siteIDTjøtta 5.656e-02 7.952e-03 2.031e+01 7.113 6.23e-07
(Intercept) ***
speciesEmpetrum nigrum ***
speciesVaccinium myrtillus ***
speciesVaccinium vitis-idaea ***
DroughtTrtExt (90)
ageClassPioneer
siteIDTjøtta ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
Overall, the model diagnostics suggest that the model fits well, with residuals that are mostly normally distributed and no major concerns with collinearity or influential outliers. However, there might be a slight issue with homoscedasticity, as indicated by the Homogeneity of Variance plot.
plot the model
# Get fitted values directly from the modelfitted_values <-predict(model_main_mean_thickness_effects , re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values #to fittedplot <-ggplot(filtered_data, aes(x = species, y = mean_thickness, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="mean_thickness",caption ="Figure 1: Shows mean_thickness of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)
interactions
# species * drought treatmentmodel_species_droughtTrt <-lmer(mean_thickness ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model_species_droughtTrt)
# Subset the data for the specified species with young and old leavessubset_data <-subset(droughtnet_data2_clean, species %in%c("Vaccinium vitis-idaea", "Empetrum nigrum"))# View the first few rows of the subsetted dataview(subset_data)
Calluna shoot type (Long and Short) only southern site had SS and LS leaves
#create a data frame for calluna lygra sitecalluna_data <- droughtnet %>%filter(species =="Calluna vulgaris"& siteID =="Lygra")view(subset_data)
# Create a data frame for Calluna vulgaris at the Lygra site and remove rows with NA valuescalluna_shoot <- droughtnet_data2_clean %>%filter(species =="Calluna vulgaris"& siteID =="Lygra") %>%drop_na()view(calluna_shoot)
SLA
#calluna shoot type * drought treatmentshoot_type_sla1 <-lmer(log(SLA) ~ calluna_shoot_type + (1| DroughNet_plotID/plant_nr), data = calluna_shoot)summary(shoot_type_sla1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(SLA) ~ calluna_shoot_type + (1 | DroughNet_plotID/plant_nr)
Data: calluna_shoot
REML criterion at convergence: -184.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.6848 -0.5395 0.1254 0.6447 2.2287
Random effects:
Groups Name Variance Std.Dev.
plant_nr:DroughNet_plotID (Intercept) 0.013066 0.11430
DroughNet_plotID (Intercept) 0.001308 0.03616
Residual 0.015593 0.12487
Number of obs: 195, groups:
plant_nr:DroughNet_plotID, 36; DroughNet_plotID, 12
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.34701 0.02537 14.71672 171.354 <2e-16 ***
calluna_shoot_typeShort -0.04379 0.01816 160.95661 -2.411 0.017 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
clln_sht_tS -0.373
---title: Analysisauthor: akformat: html: self-contained: true code-tools: true source: true fig-width: 10 fig-height: 6---```{r}#Load librarieslibrary(tidyverse)library(lme4)library(ggfortify)library(lmerTest)library(performance)library(readr)#install.packages("broom.mixed")library(broom.mixed)library(grid)```**Load data**```{r}droughtnet <-read_csv("droughtnet_data_cleaned_final.csv", col_types =cols(calluna_shoot_type =col_factor(levels =c("Long", "Short"))),guess_max =Inf)View(droughtnet)``````{r}# Specify the columns to check for NA values# List of columns to check for NA valuescolumns_to_check <-c("DroughtTrt", "DroughNet_plotID", "ageClass", "siteID", "species", "plant_height", "SLA", "mean_thickness", "dry_mass_g_original", "LDMC", "leaf_area", "wet_mass_g", "dry_mass_g")# Check for rows with any NA values in the specified columnsna_rows <-apply(droughtnet[columns_to_check], 1, function(row) any(is.na(row)))# Create a data frame showing how many NA values there are in each rowna_row_summary <- droughtnet[na_rows, columns_to_check]na_row_summary$row_na_count <-rowSums(is.na(na_row_summary))# Summarize the number of rows by the count of NAs they containna_row_summary <- na_row_summary %>%group_by(row_na_count) %>%summarise(n =n(), .groups ='drop')# Create a summary showing the co-occurrence of NA values across the specified columnsna_co_occurrence <- droughtnet %>%summarise(across(all_of(columns_to_check), ~sum(is.na(.)))) %>%pivot_longer(cols =everything(), names_to ="variable", values_to ="na_count") %>%arrange(desc(na_count))# Print the summariesprint(paste("Total rows with any NA values across checked columns:", sum(na_rows)))print(na_row_summary)print(na_co_occurrence)``````{r}# Create the new data frame with selected columnsdroughtnet_data2 <-select(droughtnet, "envelope_ID", "siteID", "species", "ageClass", "DroughtTrt", "DroughNet_plotID", "leaf_age", "calluna_shoot_type", "plant_height", "mean_thickness", "SLA", "LDMC", "plant_nr")view(droughtnet_data2)# Specify the columns to check for NA valuescolumns_to_check <-c("plant_height", "mean_thickness", "SLA", "LDMC")# Remove rows with any NA values in the specified columnsdroughtnet_data2_clean <- droughtnet_data2[!rowSums(is.na(droughtnet_data2[columns_to_check])), ]view(droughtnet_data2_clean)# Calculate the total number of NAs removed (for the specified columns only)# This calculation considers rows removed rather than individual NAstotal_NAs_removed <-nrow(droughtnet_data2) -nrow(droughtnet_data2_clean)# Calculate the total number of entries remainingtotal_entries_remaining <-nrow(droughtnet_data2_clean)# Print out the total NAs removed and total entries remainingcat("Total NAs removed:", total_NAs_removed, "\n")cat("Total entries remaining:", total_entries_remaining, "\n")```**Analysis of young leaves for the four focal species**```{r}# Specify the species of interest, including Calluna vulgarisspecies_of_interest <-c("Empetrum nigrum", "Vaccinium myrtillus", "Vaccinium vitis-idaea", "Calluna vulgaris")# Filter the data to include only the specified species where leaf_age is 'young',# and for Calluna vulgaris, additionally check that calluna_shoot_type is 'short'filtered_data <- droughtnet_data2_clean %>%filter(species %in% species_of_interest, leaf_age =="young") %>%filter(# Include Calluna vulgaris from Tjøtta regardless of calluna_shoot_type (species =="Calluna vulgaris"& siteID =="Tjøtta") |# Include Calluna vulgaris from Lygra only if calluna_shoot_type is Short (species =="Calluna vulgaris"& siteID =="Lygra"& calluna_shoot_type =="Short") |# Include other species of interest from both sites (species !="Calluna vulgaris") )view(filtered_data)```**D. plots for young leaves for the traits**```{r}filtered_data <- filtered_data %>%filter(siteID %in%c("Tjøtta", "Lygra"), DroughtTrt %in%c("Amb (0)", "Ext (90)"), ageClass %in%c("Pioneer", "Mature"))ggplot(filtered_data, aes(x = species, y = plant_height, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +# Removed leaf_age from facet_gridtheme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="Plant Height vs Treatment by Species for Young Leaves",x ="Species",y ="Plant Height") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey")) # Assigning colors to treatments``````{r}ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, vjust =0.5, hjust =1), # vertical adjustment for x axis labelsplot.caption =element_text(hjust =0.5)) +# Center the captionlabs(title ="SLA vs Treatment by Species",x ="Species",y ="Log (SLA)",fill ="Drought Treatment",caption ="Figure 2: SLA Variation in Dwarf Shrubs Under Drought Conditions Across Different Ages and Sites") +scale_fill_manual(values =c("Amb (0)"="orange", "Ext (90)"="darkgrey")) # Assigning colors to treatments``````{r}ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, vjust =0.5, hjust =1), # vertical adjustment for x axis labelsplot.caption =element_text(hjust =0.5)) +# Center the captionlabs(title ="LDMC vs Treatment by Species",x ="Species",y ="LDMC",fill ="Drought Treatment",caption ="Figure 3: LDMC Variation in dwarf Shrubs Under drought conditions across different successional phases and sites") +scale_fill_manual(values =c("Amb (0)"="darkblue", "Ext (90)"="lightblue")) # Assigning colors to treatments``````{r}filtered_data <- filtered_data %>%filter(siteID %in%c("Tjøtta", "Lygra"), DroughtTrt %in%c("Amb (0)", "Ext (90)"), ageClass %in%c("Pioneer", "Mature"))ggplot(filtered_data, aes(x = species, y = mean_thickness, fill = DroughtTrt)) +geom_boxplot() +facet_grid(ageClass ~ siteID, scales ="free") +# Removed leaf_age from facet_gridtheme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="Plant Height vs Treatment by Species for Young Leaves",x ="Species",y ="mean_thickness") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey")) # Assigning colors to treatments```**`Trait 1- LDMC- models`**```{r}# fitting a linear mixed-effects model using lmer, plant nr is the nested random effect# Ensure variables are treated as categorical if they represent categoriesfiltered_data$DroughNet_plotID <-as.factor(filtered_data$DroughNet_plotID)# Model with all the main effects for LDMCmodel_main_effects_LDMC <-lmer(LDMC ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects_LDMC)# The intercept for LDMC, representing calluna at reference level, is significantly high at 359.3762, indicating a strong baseline effect on LDMC.# Species Impact: The species Vaccinium vitis-idaea significantly reduces LDMC by 30.0715 compared to the baseline species- unique biological characteristics.# Drought Treatment: Drought treatment does not significantly affect LDMC, as indicated by a p-value of 0.160270, suggesting drought stress may not be a critical factor for LDMC variation.# Age Class Influence: Plants classified under the pioneer age class show a significant increase in LDMC (26.5932), pointing towards age or developmental stage as a determinant of LDMC.# Site Variation: Tjøtta, does not significantly influence LDMC, with a p-value of 0.152627, indicating that location-based environmental factors may not play a significant role in LDMC variationperformance::check_model(model_main_effects_LDMC, detrend =FALSE)```#The Posterior Predictive Check looks reasonable, as the observed data seem to match the model-predicted data well. #The Linearity plot suggests a linear relationship \# The Homogeneity of Variance plot does not show an obvious fan shape which would indicate heteroscedasticity, so it seems the assumption of homoscedasticity is met. \# The Influential Observations plot indicates a few points with high leverage or large residuals,-but not worrisome as they are extreme values \# The Collinearity plot suggests no issues with collinearity since all VIF values are well below the common threshold. \# The Normality of the Residuals plot indicates some deviation from normality, especially in the tails. \# The Normality of Random Effects plots show some slight deviations from normality #overall, the model fits the data well**fitting model to the ggplot**```{r}# Get fitted values directly from the modelfitted_values <-predict(model_main_effects_LDMC, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)```**interaction models**```{r}#species * DroughtTrtmodel_species_droughtTrt <-lmer(LDMC ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model_species_droughtTrt)# intercept(calluna) LDMC at 365.795 indicate a significant base level of LDMC before considering species variation or drought effects.# Species-Specific Effects:# Empetrum nigrum notably increases LDMC by 23.225 units compared to calluna demonstrating- could be a positive adaptation # Vaccinium myrtillus also shows a positive but not statistically significant effect on LDMC, suggesting a mild or variable adaptation compared to calluna# Vaccinium vitis-idaea significantly reduces LDMC by 29.776 units, highlighting its distinct physiological or ecological traits that might contribute to lower LDMC# Drought Treatment Response:# drought trt does not significantly alter LDMC on its own# Interaction with Drought:# Both Empetrum nigrum and Vaccinium myrtillus exhibit significant negative interactions with drought treatment, indicating a decrease in LDMC under drought condition- suggests these species could be sensitive to drought, affecting their LDMC.# Vaccinium vitis-idaea shows no significant interaction with drought#validate modelperformance::check_model(model_species_droughtTrt, detrend =FALSE)```**Fitting model to plot**```{r}# Get fitted values directly from the modelfitted_values <-predict(model_species_droughtTrt, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_values# Now create your plotplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}# species * age Class (h2)#this will help adress the hypothesis by examining weather the effect of drought on trait varies between age classes#significant interaction supports the hypothesis that mat and pio respond differently to droght conditions, potentially more consrvative traits in mat than pioLDMC_species_ageClass <-lmer(LDMC ~ species * ageClass + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(LDMC_species_ageClass)#validate modelperformance::check_model(LDMC_species_ageClass, detrend =FALSE)```**fitt the fitted values onto the plot**```{r}# Get fitted values directly from the modelfitted_values <-predict(LDMC_species_ageClass, re.form =NA) filtered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="lightblue", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 3: Shows overlaid points representing the fitted values of LDMC (Leaf Dry Matter Content) for each species under ambient and extreme drought treatments in two distinct successsional phase and sites.") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}# Drought Treatment * siteID (related to h3)#will help to show impact between north and south# significant interaction would indicate regional variability in how drought influences trait- hence supportinng the hypothesis, otherwise then..model1_species_siteID <-lmer(LDMC ~ species * siteID + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model1_species_siteID)#validate modelperformance::check_model(model1_species_siteID, detrend =FALSE)``````{r}# Get fitted values directly from the modelfitted_values <-predict(model1_species_siteID, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values plot <-ggplot(filtered_data, aes(x = species, y = LDMC, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="LDMC vs Treatment by Species",x ="Species", y ="LDMC",caption ="Figure 1: Shows LDMC of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}model_interactions <-lmer(LDMC ~ (species + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_interactions)#validate modelperformance::check_model(model_interactions, detrend =FALSE)```**Trait 2- plant height**```{r}#model with all the main effectsmodel_main_effects <-lmer(plant_height ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects)performance::check_model(model_main_effects, detrend =FALSE)#the model does not show any normality of residuals and also there is heterogeneity of variance#need to modify by transforming the fitted values#plant heightneeds to be log transformed to achieve homoscedasticity and linearitymodel_main_effects3 <-lmer(log(plant_height) ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects3)# Species: Empetrum nigrum, Vaccinium myrtillus, and Vaccinium vitis-idaea all show significant decreases in height compared to Calluna, indicating notable species-specific differences.# Drought Treatment: drought positively affects height, suggesting an increase under drought conditions for all species# Age Class: The pioneer age class is associated with a significant decrease in height, indicating lower levels compared to more mature.# Site: The Tjøtta site is linked to a decrease in height, highlighting the influence of location on height#validate the model- check for linearity, homoscedasticity, outliers by plotting the resudualsperformance::check_model(model_main_effects3, detrend =FALSE)```- The diagnostic plots indicate that the model has a generally good fit. The predictions match the observed data well, the residuals show linearity and homogeneity of variance, though slight deviation, and they are also normally distributed. The random effects appear to be normally distributed**fitt model to plot**```{r}# Get fitted values directly from the modelfitted_values <-predict(model_main_effects3, re.form =NA) # 'NA' to exclude random effects# Exponentiate the fitted values to transform back to the original scalefiltered_data$fitted <-exp(fitted_values)plot <-ggplot(filtered_data, aes(x = species, y = plant_height, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="Plant Height vs Treatment by Species",x ="Species", y ="Plant Height", # Updated y-axis labelcaption ="Figure 1: Shows plant height of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer and mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)```**interaction models**```{r}#models for two way interactions#species * drought trtmodel_species_drought2 <-lmer(log(plant_height) ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_species_drought2)# Compared to Calluna, the other species significantly reduce the height, with Vaccinium species showing the largest decreases.# Drought treatment on its own does not significantly affect theight, but its interaction with species (Vaccinium myrtillus and Vaccinium vitis-idaea) leads to significant increases, suggesting these species may respond positively to drought conditions in terms of thier height#plot diagnostic plotsperformance::check_model(model_species_drought2, detrend =FALSE)``````{r}# Get fitted values directly from the modelfitted_values <-predict(model_species_drought2, re.form =NA) # Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_valuesplot <-ggplot(filtered_data, aes(x = species, y =log(plant_height), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="red", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="plant height vs Treatment by Species",x ="Species", y ="LDMC") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}#species *age classmodel_height_species_ageclass <-lmer(log(plant_height) ~ species * ageClass + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_height_species_ageclass)#plot diagnostic plotsperformance::check_model(model_height_species_ageclass, detrend =FALSE)``````{r}# Get fitted values directly from the modelfitted_values <-predict(model_height_species_ageclass, re.form =NA) # Add the fitted values to the filtered_data dataframefiltered_data$.fitted <- fitted_values# Now create your plotplot <-ggplot(filtered_data, aes(x = species, y =log(plant_height), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = .fitted), position =position_dodge(width =0.75), color ="red", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1)) +labs(title ="plant height vs Treatment by Species",x ="Species", y ="plant height") +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}#drought * siteIDmodel_species_siteID <-lmer(log(plant_height) ~ species * siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_species_siteID)#plot diagnostic plotsperformance::check_model(model_species_siteID, detrend =FALSE)``````{r}model_three_interactions <-lmer(log(plant_height) ~ (species + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_three_interactions)#validate modelperformance::check_model(model_three_interactions, detrend =FALSE)```**Trait 3- SLA**```{r}# model for sla with all the fixed effectsmodel_main_effects <-lmer(SLA ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects)performance::check_model(model_main_effects, detrend =FALSE)##the model does not fit the data well#need to transform the observed values to achive normality and homoogeneity of variance#log transform SLA to achive normality and homodeneity of variancesmodel_main_effects4 <-lmer(log(SLA) ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_main_effects4)# Calluna exhibits a significant baseline SLA value of 4.327885, serving as the reference for other species comparisons.# Species Effects:# Empetrum nigrum Significantly increases SLA by 0.519248 compared to Calluna.# Vaccinium myrtillus Shows a substantial increase in SLA by 1.195279, the largest among the species.# Vaccinium vitis-idaea also significantly increases SLA by 0.513608, similar to Empetrum nigrum.# All species effects are highly significant, indicating distinct SLA characteristics compared to Calluna.# Drought treatment effect on SLA is negligible (-0.001773) and not statistically significant, suggesting drought has minimal/no impact on SLA across the specie# Age Class (Pioneer)indicates a significant decrease in SLA (-0.083333), suggesting pioneer plants have a lower SLA compared to more mature stages.# Site (Tjøtta): Associated with a significant decrease in SLA (-0.152224), indicating environmental or locational factors at Tjøtta contribute to lower SLA values#validate the model- check for linearity, homoscedasticity, outliers by plotting the resudualsperformance::check_model(model_main_effects4, detrend =FALSE)```Overall, the model shows a generally good fit with some concerns:- The posterior predictive check indicates a good fit. The linearity assumption seems to be met. There is a potential issue with homoscedasticity, as the Homogeneity of Variance plot suggests that variance slightly increases with fitted values. Influential observations do not appear to be a problem. The residuals and random effects are approximately normally distributed**fit model onto the plot**```{r}# Get fitted values directly from the modelfitted_values <-predict(model_main_effects4, re.form =NA) # Exponentiate the fitted values to transform back to the original scalefiltered_data$fitted <-exp(fitted_values)# Now create your plot with a caption at the bottomplot <-ggplot(filtered_data, aes(x = species, y = SLA, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +# Use exponentiated fitted valuesfacet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA", caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer and mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)```**Interaction models for SLA**```{r}#species * droughtTrtsla_species_droughtTrt <-lmer(log(SLA) ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(sla_species_droughtTrt)performance::check_model(sla_species_droughtTrt, detrend =FALSE)```In summary, the model seems to have a good fit overall. high VIF for the 'species:DroughtTrt' interaction term is acceptable, and slight deviations from normality in the residuals. As for the residuals, the minor deviations from normality .```{r}# Get fitted values directly from the modelfitted_values <-predict(sla_species_droughtTrt, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))print(plot)``````{r}#species * ageclasssla_species_ageClass <-lmer(log(SLA) ~ species * ageClass + (1| DroughNet_plotID), data = filtered_data)summary(sla_species_ageClass)performance::check_model(sla_species_ageClass, detrend =FALSE)``````{r}# Get fitted values directly from the modelfitted_values <-predict(sla_species_ageClass, re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <-exp(fitted_values) # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y = SLA, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)``````{r}# species * siteIDsla_species_siteID <-lmer(log(SLA) ~ species * siteID + (1| DroughNet_plotID), data = filtered_data)summary(sla_species_siteID)performance::check_model(sla_species_siteID, detrend =FALSE)``````{r}# Get fitted values directly from the modelfitted_values <-predict(sla_species_siteID, re.form =NA) # 'NA' to exclude random effectsfiltered_data$fitted <- fitted_values # Corrected .fitted to fittedplot <-ggplot(filtered_data, aes(x = species, y =log(SLA), fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="SLA",caption ="Figure 1: Shows SLA of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)``````{r}model_three_interactions <-lmer(log(SLA) ~ (species + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_three_interactions)#validate modelperformance::check_model(model_three_interactions, detrend =FALSE)```**Trait 4- mean_thickness**```{r}#general model#model with all the main effectsmodel_main_mean_thickness_effects <-lmer(mean_thickness ~ species + DroughtTrt + ageClass + siteID + (1| DroughNet_plotID/plant_nr),data = filtered_data) summary(model_main_mean_thickness_effects )performance::check_model(model_main_mean_thickness_effects , detrend =FALSE)```Overall, the model diagnostics suggest that the model fits well, with residuals that are mostly normally distributed and no major concerns with collinearity or influential outliers. However, there might be a slight issue with homoscedasticity, as indicated by the Homogeneity of Variance plot.**plot the model**```{r}# Get fitted values directly from the modelfitted_values <-predict(model_main_mean_thickness_effects , re.form =NA) # 'NA' to exclude random effects# Add the fitted values to the filtered_data dataframefiltered_data$fitted <- fitted_values #to fittedplot <-ggplot(filtered_data, aes(x = species, y = mean_thickness, fill = DroughtTrt)) +geom_boxplot() +geom_point(aes(y = fitted), position =position_dodge(width =0.75), color ="brown", size =1.5) +facet_grid(ageClass ~ siteID, scales ="free") +theme_bw() +theme(axis.text.x =element_text(angle =90, hjust =1),plot.margin =unit(c(1, 1, 2, 1), "lines"), plot.caption =element_text(hjust =0.5) # Center align the caption ) +labs(title ="SLA vs Treatment by Species",x ="Species", y ="mean_thickness",caption ="Figure 1: Shows mean_thickness of the four species under ambient and extreme drought conditions in two distinct successional phases (pioneer nad Mature) and two distinct sites (Tjøtta in the north and Lygra in the south)." ) +scale_fill_manual(values =c("Amb (0)"="blue", "Ext (90)"="grey"))# Print the plotprint(plot)```**interactions**```{r}# species * drought treatmentmodel_species_droughtTrt <-lmer(mean_thickness ~ species * DroughtTrt + (1| DroughNet_plotID/plant_nr), data = filtered_data)summary(model_species_droughtTrt)performance::check_model(model_species_droughtTrt, detrend =FALSE)``````{r}#species * ageclassmodel_species_ageClass <-lmer(mean_thickness ~ species * ageClass + (1| DroughNet_plotID), data = filtered_data)summary(model_species_ageClass)performance::check_model(model_species_ageClass, detrend =FALSE)``````{r}#species * ageclassmodel_species_ageClass <-lmer(mean_thickness ~ species * ageClass + (1| DroughNet_plotID), data = filtered_data)summary(model_species_ageClass)performance::check_model(model_species_ageClass, detrend =FALSE)``````{r}#species * siteIDmodel_species_siteID <-lmer(mean_thickness ~ species * siteID + (1| DroughNet_plotID), data = filtered_data)summary(model_species_siteID)#validate the modelperformance::check_model(model_species_siteID, detrend =FALSE)``````{r}model_three_interactions <-lmer(mean_thickness ~ (species + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr),data = filtered_data)summary(model_three_interactions)#validate modelperformance::check_model(model_three_interactions, detrend =FALSE)```**Analysis for leaf age (young and old leaves)**```{r}# Subset the data for the specified species with young and old leavessubset_data <-subset(droughtnet_data2_clean, species %in%c("Vaccinium vitis-idaea", "Empetrum nigrum"))# View the first few rows of the subsetted dataview(subset_data)```**SLA**```{r}#1 slasla_leafAge <-lmer(log(SLA) ~ leaf_age + (1| DroughNet_plotID/plant_nr), data = subset_data)summary(sla_leafAge)performance::check_model(sla_leafAge, detrend =FALSE)``````{r}sla_two_way_interactions <-lmer(log(SLA) ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^2+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(sla_two_way_interactions)performance::check_model(sla_two_way_interactions, detrend =FALSE)``````{r}sla_three_way_interactions <-lmer(log(SLA) ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(sla_three_way_interactions)performance::check_model(sla_three_way_interactions, detrend =FALSE)```**LDMC**```{r}#2 ldmcldmc_leafAge <-lmer(LDMC ~ leaf_age + (1| DroughNet_plotID/plant_nr), data = subset_data)summary(ldmc_leafAge)performance::check_model(ldmc_leafAge, detrend =FALSE)``````{r}ldmc_two_way_interactions <-lmer(LDMC ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^2+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(ldmc_two_way_interactions)performance::check_model(ldmc_two_way_interactions, detrend =FALSE)``````{r}ldmc_three_way_interactions <-lmer(log(SLA) ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(ldmc_three_way_interactions)performance::check_model(ldmc_three_way_interactions, detrend =FALSE)```**mean thickness**```{r}species_leafAge <-lmer(mean_thickness ~ species * leaf_age + (1| DroughNet_plotID/plant_nr), data = subset_data)summary(species_leafAge)performance::check_model(species_leafAge, detrend =FALSE)``````{r}thickness_two_way_interactions <-lmer(mean_thickness ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^2+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(thickness_two_way_interactions)performance::check_model(thickness_two_way_interactions, detrend =FALSE)``````{r}thickness_three_way_interactions <-lmer(mean_thickness ~ (species + leaf_age + DroughtTrt + ageClass + siteID)^3+ (1| DroughNet_plotID/plant_nr), data = subset_data)summary(thickness_three_way_interactions)performance::check_model(thickness_three_way_interactions, detrend =FALSE)```**Calluna shoot type (Long and Short) only southern site had SS and LS leaves**```{r}#create a data frame for calluna lygra sitecalluna_data <- droughtnet %>%filter(species =="Calluna vulgaris"& siteID =="Lygra")view(subset_data)``````{r}# Create a data frame for Calluna vulgaris at the Lygra site and remove rows with NA valuescalluna_shoot <- droughtnet_data2_clean %>%filter(species =="Calluna vulgaris"& siteID =="Lygra") %>%drop_na()view(calluna_shoot)```**SLA**```{r}#calluna shoot type * drought treatmentshoot_type_sla1 <-lmer(log(SLA) ~ calluna_shoot_type + (1| DroughNet_plotID/plant_nr), data = calluna_shoot)summary(shoot_type_sla1)performance::check_model(shoot_type_sla1, detrend =FALSE)``````{r}#calluna shoot typeshoot_type_sla2 <-lmer(log(SLA) ~ (calluna_shoot_type + DroughtTrt + ageClass)^2+ (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type_sla2)performance::check_model(shoot_type_sla2, detrend =FALSE)``````{r}#calluna shoot typeshoot_type_sla3 <-lmer(log(SLA) ~ (calluna_shoot_type * DroughtTrt * ageClass) + (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type_sla3)performance::check_model(shoot_type_sla3, detrend =FALSE)```**LMDC**```{r}shoot_type1 <-lmer(LDMC ~ calluna_shoot_type + (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type1)performance::check_model(shoot_type1, detrend =FALSE)``````{r}shoot_type_lmdc2 <-lmer(LDMC ~ (calluna_shoot_type + DroughtTrt + ageClass)^2+ (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type_lmdc2)performance::check_model(shoot_type_lmdc2, detrend =FALSE)``````{r}#calluna shoot typeshoot_type_lmdc3 <-lmer(LDMC ~ (calluna_shoot_type * DroughtTrt * ageClass) + (1| DroughNet_plotID/plant_nr), data = calluna_shoot)summary(shoot_type_lmdc3)performance::check_model(shoot_type_lmdc3, detrend =FALSE)```**Thickness**```{r}shoot_type_thickness1 <-lmer(mean_thickness ~ calluna_shoot_type + (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type_thickness1)performance::check_model(shoot_type_thickness1, detrend =FALSE)``````{r}shoot_type_thickness2 <-lmer(mean_thickness ~ (calluna_shoot_type + DroughtTrt + ageClass)^2+ (1| DroughNet_plotID), data = calluna_shoot)summary(shoot_type_thickness2)performance::check_model(shoot_type_thickness2, detrend =FALSE)``````{r}shoot_type_thickness3 <-lmer(mean_thickness ~ (calluna_shoot_type * DroughtTrt * ageClass) + (1| DroughNet_plotID/plant_nr), data = calluna_shoot)summary(shoot_type_thickness3)performance::check_model(shoot_type_thickness3, detrend =FALSE)```
